Robust circadian clocks from coupled protein modification and 
transcription-translation cycles 
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The cyanobacterium Synechococcus elongatus uses both a protein phosphorylation cycle and a 
transcription-translation cycle to generate circadian rhythms that are highly robust against bio- 
chemical noise. We use stochastic simulations to analyze how these cycles interact to generate 
stable rhythms in growing, dividing cells. We find that a protein phosphorylation cycle by itself is 
robust when protein turnover is low. For high decay or dilution rates (and co mpensating synthe- 
sis rate), however, the phosphorylation-based oscillator loses its integrity. Circadian rhythms thus 
cannot be generated with a phosphorylation cycle alone when the growth rate, and consequently 
the rate of protein dilution, is high enough; in practice, a purely post-translational clock ceases to 
function well when the cell doubling time drops below the 24 hour clock period. At higher growth 
rates, a transcription-translation cycle becomes essential for generating robust circadian rhythms. 
Interestingly, while a transcription-translation cycle is necessary to sustain a phosphorylation cy- 
cle at high growth rates, a phosphorylation cycle can dramatically enhance the robustness of a 
transcription-translation cycle at lower protein decay or dilution rates. Our analysis thus predicts 
that both cycles are required to generate robust circadian rhythms over the full range of growth 
conditions. 
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I. INTRODUCTION 

Many organisms use circadian clocks to anticipate 
changes between day and night It had long 

been believed that these clocks are driven primarily by 
transcription-translation cycles (TTCs) built on nega- 
tive feedback. However, while some circadian clocks can 
maintain robust rhythms for years in the absence of any 
daily cue recent experiments have vividly demon- 
strated that gene expression is often highly stochastic 
. This raises the question of how these clocks can be 
robust against biochemical noise. In multicellular organ- 
isms, the robustness might be explained by intercellular 
interactions 0, H| j but it is now known that even unicel- 
lular organisms can have very stable circadian rhythms. 
The clock of the cyanbacterium Synechococcus elongatus, 
for example, has a correlation time of several months Q , 
even though the clocks of the different cells in a popu- 
lation hardly interact with one another |6, 9]. Interest- 
ingly, it has recently been discovered that the S. elon- 
gatus clock also includes a protein phosphorylation cycle 
(PPC) that can run independently of the TTC 0, ij. It 
has been suggested that this protein modification oscil- 
lator [8] . possibly in combination with the transcription- 
translation oscillator [l5j], is responsible for the clock's 
striking noise resistance. Here, we use mathematical 
modeling to study how these two clocks interact in grow- 
ing, dividing cells. We find that the PPC alone is very 
stable when cell growth is negligible, but it begins to fail 
as the growth rate increases, and for doubling times typ- 
ical of S. elongatus it cannot function without help from 
a TTC. We then use stochastic simulations to show that 
a PPC can, however, dramatically enhance the robust- 
ness of a TTC, especially when cells are growing slowly. 
The two mechanisms thus perform best in complemen- 



tary situations, and it is likely that both oscillators are 
necessary to maintain even minimal clock stability across 
the full range of conditions and growth rates that bacte- 
ria encounter in the wild. Importantly, the TTC and the 
PPC in S. elongatus are much more tightly intertwined 
than conventional coupled phase oscillators; as a result, 
the combination of the two far outperforms not just each 
of its two components individually, but also a hypotheti- 
cal system in which the two parts are coupled in normal 
textbook fashion [Icl |. 

In S. elongatus, the central components of the clock 
are the three genes kaiA, kaiB, and kaiC [llj]. Under 
continuous light conditions, the levels of mRNA from the 
kaiBC operon and of the protein KaiC oscillate in a cir- 
cadian fashion with a 6 hour lag between their peaks [l2j . 
Overexpression of phosphorylated KaiC abolishes kaiBC 
expression [TH, [l4| , while a transient increase in KaiC re- 
sets the phase of the oscillator [lS^ . These observations 
led to the proposal that the Kai system is a transcription- 
translation oscillator, with KaiC negatively regulating its 
own transcription. In 2005, however, Kondo and cowork- 
ers showed that KaiC is phosphorylated in a cyclical man- 
ner with a period of 24 hours, even when kaiBC transcrip- 
tion is inhibited Q • Still more remarkably, the rhythmic 
phosphorylation of KaiC could be reconstituted in the 
test tube in the presence of only KaiA, KaiB, and ATP 
Q. This raised the possibility that the principal pace- 
maker of the clock is not a transcription-translation cycle, 
but a protein phosphorylation cycle Yet, in 2008, the 
same group showed that circadian oscillations of gene ex- 
pression persist even when KaiC is always held in a highly 
phosphorylated state [l5[ . They thus concluded that the 
clock is driven by both a TTC and a PPC, and suggested 
that the interactions between the two oscillators may en- 
hance the robustness of the clock 15 1. 
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The PPC has been characterized experimentally in 
considerable detail. It is known, for exam ple, that KaiC 
forms a homo-hexamer [l5[ , KaiA a dimer |l5| , and KaiB 
a dimer or a tetramer [jj} ■ KaiC has two phospho- 

rylation sites per protein monomer, which are phospho- 
rylated and dephosphorylated in a definite sequence as a 
result of KaiC's autokinase and autophosphatase activ- 
ity 0, HI • KaiA stimulates KaiC phosphorylation [20l . l2~fl | , 
while KaiB negates the effect of KaiA [2(^11. Thanks 
to the wealth of available experimental data, the PPC 
has proven an unusually fruitful system for mathemati- 
cal modeling [l|, 0, EE MM, ■ 

The TTC appears to encompass several different regu- 
latory pathways and is much less well understood. Multi- 
ple lines of evidence indicate that circadian promoter ac- 
tivity can be achieved without any specific cis-regulatory 
element [33l - l35l j . suggesting that the clock may regu- 
late transcription at least in part through changes in 
chromosome compaction and superhelicity [H|, EE 
Several proteins important for transcriptional regula- 
tion of the kaiBC operon have also been identified. In 
particular, LabA and the sensor histidine kinase CikA 
have been implicated in negative feedback of phospho- 
rylated KaiC on its own transcription during the sub- 
jective night [U El EI El El, while the SasA histi " 
dine kinase, acting through the RpaA response regulator, 
seems to mediate positive feedback during the subjec- 
tive day [H SU . Epistasis analysis suggests that LabA 
may likewise act through RpaA ;38f. The picture that 
emerges is thus that during the phosphorylation phase of 
the PPC, KaiC activates kaiBC expression, while in the 
dephosphorylation phase, it represses kaiBC expression 

M (Fig- P- 

In this study, we perform stochastic simulations to in- 
vestigate the roles of the PPC and the TTC. We first 
examine a situation in which the total number of each 
Kai protein is constant — they are neither produced nor 
destroyed — and only the PPC is operative. We show that 
in this case the PPC is highly robust against noise aris- 
ing from the intrinsic stochasticity of chemical reactions. 
Even for reaction volumes smaller than the typical vol- 
ume of a cyanobacterium, the correlation time is longer 
than that observed experimentally [9]. Living cells, how- 
ever, constantly grow and divide, and proteins must thus 
be synthesized to balance dilution. In fact, dilution can 
be thought of as introducing an effective protein degra- 
dation rate set by the cell doubling time. We therefore 
next study a PPC in which the Kai proteins are produced 
and degraded with rates that are constant in time. The 
simulations reveal that protein synthesis and decay dra- 
matically reduce the viability of the PPC; we predict that 
for a cell doubling time of 24 hours and a bacterial volume 
of 1 /im 3 , the PPC dephases in roughly 10 days, much 
faster than real S. elongatus [9J]. The same loss of oscil- 
lations is seen even in the limit of infinite volume, indi- 
cating it is not due solely to noise. Instead, the constant 
synthesis of proteins, which we assume are all initially 
created in the same phosphorylation state, necessarily in- 




FIG. 1: A cartoon of how the TTC might interact with the 
PPC [H. Active RpA activates kaiBC expression. During 
the phosphorylation phase, KaiC activates RpaA, while dur- 
ing the dephosphorylation phase, KaiC represses RpaA. 



jects KaiC hexamers with the "wrong" phosphorylation 
levels at certain phases of the cycle [3, HI; if these appear 
fast enough, they can destroy the oscillation. One role 
of the TTC is thus to introduce proteins only when the 
phosphorylation state of the freshly made KaiC matches 
that of the PPC. Our simulations of a combined PPC and 
TTC reveal that a TTC can indeed greatly enhance the 
robustness of the PPC, yielding correlation times consis- 
tent with those measured experimentally [9] . Finally, we 
consider whether the PPC is needed at all, or whether one 
could build an equally good circadian clock using only a 
TTC. We find that it is possible to construct a TTC with 
a period of 24 hours and the observed correlation time 
of a few months [9J] . However, this comes at the expense 
of very high protein synthesis and decay rates, which im- 
pose an extra energetic burden on the cell. Our results 
thus suggest that a PPC allows for a more robust oscilla- 
tor at a lower cost. Although our models are simplified, 
we argue in the Discussion section that our qualitative 
results are unavoidable consequences of the interaction 
between a circadian clock and cell growth and so should 
hold far more generally. 



II. RESULTS 

A. A protein phosphorylation cycle with constant 
protein concentrations is highly robust 

In recent years, the PPC has been mathematically 
modelled with considerable success [H, 0, EE EHM & 
El| (for a review, see [11). All models endow KaiC with 
an intrinsic ability to cyclically phosphorylate and de- 
phosphorylate itself, but they differ in how the cycles of 
the individual KaiC proteins are synchronized. In this 
manuscript, we adopt the model developed by us [l|, 
which is one of several based on a mechanism that we 
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called "differential affinity" 0, H|: KaiA stim- 
ulates KaiC phosphorylation, but the limited supply of 
KaiA dimers binds preferentially to those KaiC molecules 
that are falling behind in the cycle, allowing them to 
catch up. Specifically, in our model each KaiC hexamer 
can switch between an active conformational state C,, 
where the number i of phosphorylated monomers tends 
to increase, and an inactive state Gi, where i tends to 
decrease (Fig. [1]); KaiA stimulates phosphorylation of 
active KaiC, but is sequestered by complexes containing 
KaiB and inactive KaiC. KaiC in the inactive state can 
thus delay the progress of fully dephosphorylated hexam- 
ers that have already switched back to the active state 
and are ready to be phosphorylated again. With A and B 
denoting, respectively, a KaiA dimer and a KaiB dimer, 
the model becomes: 



G, f± Q, Q + A «=* AC, 4 Ci+i + A 

h k Ab 



Gi + B «^ BCi, BCi + B «^ B2Ci 

fe? b 2fc? b 



(1) 



(2) 



xk At k At 

B X G,+A <± AB x Ci,AB 2 Ci+A <=* A 2 B 2 Q (3) 

k Ab 2k Ab 

k k 

Cj Gi+i,Cj ^ Ci+i (4) 

&dps k dps 

h h 

B^C^ B^C^+i, AyB^Ci AyB^Ci+i . (5) 

&dps ^dps 

This model reproduces the phosphorylation behaviour 
of KaiC in vitro not only when all Kai proteins are 
present, but also when KaiA and/or KaiB are absent 0. 
It moreover correctly predicted the experimentally ob- 
served disappearance of oscillations when the KaiA con- 
centration is raised 0,0], a success that strongly supports 
the idea that KaiA sequestration is the primary driver of 
synchronization. Our model does not feature monomer 
exchange between KaiC hexamers, an alternative means 
of synchronisation [25j | that has been observed in experi- 
ments [5J, |29[ ; we and others find that monomer exchange 
is not critical for stable oscillations 

0,0, SI- Inthe Sup- 
porting Information (ST), we show that similar results are 

obtained with a model that focuses on the phosphoryla- 
tion cycle of individual KaiC monomers Q rather than 
of KaiC hexamers and that includes a positive feedback 
loop. 

We quantify our model's robustness to chemical noise 
by performing kinetic Monte Carlo simulations of the 
chemical master equation [J]. In our simulations, we 
vary the reaction volume, but keep the concentrations of 
the Kai proteins constant at levels comparable to those 
used in the in vitro experiments 0]. Fig. shows 
a time trace of the fraction p(t) of phosphorylated KaiC 
monomers at a volume of order that of a cyanobacterium, 
while Fig. 03 gives, as a function of volume, the corre- 
lation number of cycles Ui/2, defined as the number of 
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FIG. 2: The KaiC phosphorylation oscillator is highly stable 
when proteins are not synthesized or degraded. (A) Time 
trace of the phosphorylation ratio pit), i.e. the fraction of 
monomers that are phosphorylated, for volumes V = 1/xm 3 
and V = 0.1/^m 3 . (B) Correlation number of cycles, ni/2, as a 
function of the volume. The protein concentrations are those 
used in the in vitro experiments [A]t = 0.58/iM; [B]t = 

1.75/iM; [C] T = 0.58^M. For other parameters, see Table SI 
of SI. 



cycles after which the standard deviation in the phase of 
the oscillation is half a day ; in the SI we discuss our 
error estimates, which suggest that our computed n\j2 
values are typically accurate to within 15-20%. One is- 
sue that arises in comparing the simulation results to 
the measured in vivo clock robustness is that the Kai 
proteins appear to be present in living cells in a ratio 
at which the in vitro system would not oscillate. Ki- 
tayama et al. found that there are approximately 10,000 
KaiC monomers but only 250-500 KaiA monomers per S. 
elongatus cell [22|], corresponding to a roughly 6:1 KaiC 
hexamer to KaiA dimer ratio instead of the standard 
1:1 ratio in vitro 0, 0|. It has been suggested that this 
discrepancy may indicate that the clock reactions are in 
fact confined to a subdomain of the whole cell from which 
some of the KaiB and KaiC molecules are excluded [22j ; 
here, we adopt this hypothesis and assume that the Kai 
proteins are found in the physiologically relevant reac- 
tion volume in proportions comparable to those used in 
the in vitro experiments. If we take this volume to be 
~ lyum 3 , comparable to the size of the entire cell, then 
Fig. [2)3 shows that ru /% « 200, on the high side of the 
measured 166 ± 100 [9(. Even if we consider a volume of 
only 0.5/im 3 — small enough that the measured number 
of KaiA molecules is adequate to give the in vitro KaiA 
dimer concentration of 0.58/iM — we find that ni« ~ 102, 
still within the experimental bounds (in contrast to the 
predictions of some alternative models 0,113; see SI). 
Our model thus predicts that the PPC is very resistant 
to noise arising from the intrinsically stochastic nature of 
chemical reactions. 



B. A phosphorylation cycle with constant protein 
synthesis and degradation rates is not stable 

Fig. [2] shows that the phosphorylation cycle is highly 
robust when the total concentrations of the Kai proteins 
are strictly constant. But in vivo proteins are constantly 
being synthesized and degraded. To study how this af- 
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FIG. 3: A system in which the Kai proteins are produced 
and degraded with rates that are constant in time cannot 
sustain a stable phosphorylation cycle. (A) Time trace of 
the phosphorylation ratio p(t) for three different degradation 
rates n and V — 1/im 3 . (B) Correlation number of cycles, 
ni/2, as a function of the volume and degradation rate. The 
white line denotes the bifurcation line fi = 0.0621h _1 , where 
the system undergoes a supercritical Hopf bifurcation in the 
deterministic limit; the three colored dots correspond to the 
three curves in panel A. All proteins are degraded with the 
same rate; KaiA, KaiB and KaiC are produced with rates 
such that the average total concentrations equal those used 
in the in vitro experiments (see Fig. 



fects the PPC, we consider a model in which the Kai pro- 
teins are produced and degraded in a stochastic (mcmo- 
ryless) fashion with rates that are constant in time, with 
the effects of active degradation |2| and of passive dilu- 
tion lumped into a single first-order decay rate \x (see 
SI). 

Fig. [3] shows the performance of the phosphorylation 
cycle as a function of degradation rate and cell volume; 
in all cases, the synthesis rates are adjusted so that the 
mean concentrations are constant and equal to those 
used in the previous section. The oscillator's robust- 
ness clearly decreases dramatically with increasing pro- 
tein synthesis and decay rate. For a volume comparable 
to that of a cyanobacterium and a degradation rate of 
0.03hr~ , the correlation time is less than 20 days, much 
lower than that observed in vivo *9i]. This degradation 
rate is precisely the effective rate arising from protein di- 
lution with a cell doubling time of 24 hours. It is known, 
however, that KaiC is also degraded actively at a rate as 
high as O.lhr -1 Q, leading to still worse stability. 

The disappearance of the oscillations for higher protein 
synthesis and decay rates can be understood by noting 
that fresh KaiC hexamers are made in a fixed phosphory- 
lation state, which then has to catch up with those of the 
proteins that are already in the cycle 0, |42[ . When the 
degradation rate is high, the new proteins are likely to 
be degraded before the PPC can synchronise their phos- 
phorylation levels; indeed, in the limit that the protein 
synthesis and decay rates go to infinity, p{t) becomes con- 
stant in time and equal to the phosphorylation level of 
freshly made KaiC proteins. This is not a purely stochas- 
tic effect; the dashed bifurcation line of Fig. [3j3 shows 
that even in a deterministic model the oscillations disap- 
pear when the synthesis and decay rates become too big 
(see SI). 



C. A protein phosphorylation cycle with a 
transcription-translation cycle is very stable 



To sustain a phosphorylation cycle, KaiC has to be 
made in an oscillatory fashion: newly synthesized KaiC 
proteins have to be injected into the phosphorylation 
cycle at the moment that the phosphorylation states of 
the other proteins that are already in the cycle match 
the modification state of the fresh KaiC proteins [l[. 
This is the principal role of the transcription-translation 
cycle. Here, we present a model for how such a cycle 
might interact with the protein phosphorylation cycle 
(see Fig. [T]). The model is inspired by that of Kondo et 
al. [39j and contains the following key ingredients: 

1. RpaA activates kaiBC expression 

Deletion of rpaA reduces the expression of clock- 
controlled genes, including that of kaiBC [H, HH, l4l| . 
Since neither promoters nor transcription or 
chromosome-compaction factors have been identified 
that interact with RpaA [4l| . we make the phenomeno- 
logical ass ump tion that RpaA directly activates kaiBC 
expression [39[. 

2. RpaA is activated by KaiC when KaiC is in the active 
state 

RpaA is activated via phosphorylation by the histidine 
kinase SasA, whose activity is in turn stimulated by KaiC 
36j [ill: inactivation of SasA reduces kaiBC expression 
38l I39L |41|. Moreover, SasA- RpaA phosphorylation oc- 
curs 4-8 hours before the peak of KaiC phosphorylation 
41]. This suggests that partially phosphorylated KaiC 
that is on the active branch activates RpaA through 
SasA [39] . Since SasA phosphorylation, occurring on 
time scales of minutes 36], is much faster than KaiC 
phosphorylation, occurring on time scales of hours, we 
assume that the SasA dynamics can be integrated out. 

3. RpaA is inactivated by KaiC when KaiC is in the 
inactive state 

Inactivation of LabA [38[ or CikA [39[ increases kaiBC 
expression, with inactivation of both having a still 
stronger effect (39|. SasA inactivation can compen- 
sate for both LabA [38| and CikA inactivation [39j. 
but RpaA inactvation cannot compensate for LabA 
Taken together, these results suggest 
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inactivation 

that SasA, LabA, and CikA control kaiBC expression 
through different pathways, with at least the SasA and 
LabA pathways converging on RpaA (39|. Because 
phosphorylation of KaiC is critical for negative feedback 
on kaiBC expression 14 1. LabA and CikA appear to 
act downstream of phosphorylated KaiC (39|. Since the 
mechanisms by which LabA and/or CikA repress RpaA 
activation are unknown, we make the phenomenological 
assumption that KaiC on the inactive branch deactivates 
RpaA. 

4- KaiC is injected into the system as fully phosphory- 
lated hexamers 

Imai et al. reported that newly synthesized KaiC is 
phosphorylated in vivo within 30 min Q, much faster 
than phosphorylation of KaiC hexamers in vitro, which 
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takes about 6 hours [7|; we thus assume that newly 
synthesized KaiC is injected into the system as fully 
phosphorylated KaiC hexamers. 

5. KaiA and RpaA are synthesized at constant rates 
The mRNA levels of kaiA and rpaA exhibit much 
weaker oscillations than those of kaiB and kaiC [H|]. We 
therefore assume that kaiA and rpaA are expressed at 
constant rates. 

6. The phosphorylation cycle in vivo is similar to that 
in vitro. 

Fig. [1] shows a cartoon of this model, which is 
described by the reactions of Eqs. HHS] for the PPC 
together with the following reactions for the TTC and 
the coupling between them: 




Time [h] Volume V [urn 3 ] 



FIG. 4: The combination of a TTC and a PPC can gener- 
ate stable circadian rhythms. (A) Time trace of the phos- 
phorylation ratio p(t) and total KaiC concentration [C]t for 
V — l/im s and /i = 0.03hr _1 (solid lines) and O.lhr -1 (dashed 
lines). (B) Correlation number of cycles, «i/2, for p(t) as a 
function of the volume and degradation rate fj,. The protein 
synthesis rates are varied with the degradation rates such that 
the average protein concentrations equal those used in vitro 
(see Fig. d]). 



R + X*^R + X, R+X^R+X 



/3 c [R] 4 /(^ 4 +[R] 4 ) 



C 6 + 3B 



R, R, A, B, ACi, B^Ci, A^B^C^ 



(0) 
(7) 

(8) 
(9) 



Eq. [5] models activation of inactive RpaA, R, by X £ 
{AC2,...ACs} and inactivation of active RpaA, R, by 
X e {AyB x C5 . . . A y B x C2}- The model is robust to the 
precise choice of phosphoforms that activate and repress 
RpaA (see ST). Because KaiA enhances kaiBC expres- 
sion [2(1, we assume that active KaiC bound to KaiA 
activates RpaA; a model in which both Cj and AQ acti- 
vate RpaA also generates oscillations, albeit less robustly 
(see ST). Eq. [7] models the activation of kaiBC expres- 
sion by RpaA, using a Hill function with coefficient 4. We 
assume a normally distributed delay, denoted by the dou- 
ble arrow, with mean r = 5 hr and standard deviation 
<j t = 0.5 hr, between the activation of kaiBC transcrip- 
tion and the appearance of KaiB and KaiC protein. The 
length of the delay is dictated by the requirement that 
new KaiC protein be produced when its phosphorylation 
state matches that of the PPC. Eq. [S]models constitutive 
expression of kaiA and rpaA. Eq. [^describes degradation 
of all species with the same rate constant /i; we ignore 
rhythmic KaiC degradation Q, which is not essential to 
produce a robust clock. In the SI we discuss a more de- 
tailed model that includes cell growth and binomial par- 
titioning upon cell division, which appreciably increases 
noise without changing qualitative trends. Although one 
can think of our model in loose terms as consisting of cou- 
pled transcriptional and post-translational oscillators, it 
cannot be mathematically decomposed into two separate 
oscillatory systems, each with its own variables, nor can 
the strength of the coupling between the two cycles be 
independently tuned. It is thus formally quite different 
from textbook models of coupled oscillators 10]. 

Fig. S3 shows the robustness of this PPC-TTC model 
as a function of cell volume and protein degradation rate 
fi; when \x is varied, all three protein synthesis rates j3 a 



are adjusted so that the average protein concentrations 
arc unchanged and comparable to those of the models 
considered above. As expected, n\/2 decreases with de- 
creasing cell volume. Its dependence on the degradation 
rate, however, is markedly different from that seen with 
constant KaiC synthesis (Fig. 03): A PPC sustained by 
a TTC becomes more robust with increasing decay rate. 
If we assume that proteins are lost only through dilu- 
tion, then for a bacterial volume of 1 fim 3 and a cell 
doubling time of 24 hours (corresponding to a decay 
rate fi = 0.03hr _1 ), the correlation time is about 200 
days, consistent with the value measured experimentally 
9]. If proteins are also degraded actively, increasing //, 
this excellent behavior improves still further; even for 
V = 0.5/im 3 , m/2 ~ 120 for [i = 0.1hr~\ This is in 
stark contrast to the stability of a PPC without a TTC, 
which we have seen is hardly stable for a protein lifetime 
of about 10 hours (Fig. EP). 



D. A protein phosphorylation cycle dramatically 
enhances the robustness of a 
transcription-translation cycle 

Figs. [3)3 and[3j3 show that a TTC can greatly improve 
the robustness of a PPC. One might thus ask whether 
the PPC is needed at all, or whether an adequate clock 
can be built with only a TTC. To address this question, 
we modify the PPC-TTC model (Eqs. EHSJ) so that it 
consists only of a TTC: 



/3 c [R] 4 /(if 4 + [R] 4 ) 



C, 0^R, C,RA0 (10) 



R^-R, C + R^C + R 



(11) 



In the simulations, we adjust the delay r and the synthe- 
sis rate f3 for each choice of the degradation rate n such 
that the oscillation period is 24 hours and the average 
KaiC concentration is comparable to that of the other 
models considered so far (see ST). We also investigated 
the simplest possible TTC model, namely one in which 
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KaiC directly represses its own synthesis; this gave very 
similar results (see SI). 

Fig. [SJ3 shows for a bacterial volume of approximately 
1 /xm 3 the robustness of this TTC-only model as a func- 
tion of the degradation rate, together with the results 
for the PPC-TTC model (Fig. gp) and for the PPC- 
only model with constant synthesis rates (Fig. |3j3). The 
TTC-only model's behavior is the opposite of that of 
the PPC-only model; the TTC-only model is most stable 
for high degradation rates, and its robustness falls dra- 
matically when fi drops below 0.2hr~ . The combined 
TTC-PPC model, however, is robust for all degradation 
rates — its correlation time interpolates between those of 
the TTC-only model for large fj, and of the PPC-only 
model for small /i. Importantly, the relative advantage 
of the combined model is greatest when 1 //j, is of order the 
oscillator period, and this is precisely the regime where 
physiological degradation and dilution rates for S. elon- 
gatus fall 2\. Further, the combined oscillator does much 
better than either oscillator alone; in contrast, when two 
conventional phase oscillators with comparable noise lev- 
els are coupled, one expects only about a factor of two 
gain in n 1 / 2 [10|]. 

Dilution puts a lower bound on the degradation rate, 
which means that a stable oscillator cannot be based 
on a PPC only, especially when the growth rate of the 
bacterium is large. The degradation rate can, on the 
other hand, be increased by active degradation. For high 
degradation rates, i.e. fx > 0.2hr _ , an oscillator based 
only on transcriptional feedback can be sufficiently ro- 
bust (Fig. [SJ3). However, to balance these high decay 
rates, the protein synthesis rates have to be correspond- 
ingly large, which is energetically costly [13]. Combin- 
ing a PPC with a TTC makes it possible to dramatically 
improve the robustness while keeping the protein synthe- 
sis rates the same. While a p hosphorylation cycle does 
not come entirely for free [48|], this suggests that a PPC 
combined with a TTC gives the best performance-to-cost 
ratio. 

The question remains why a PPC so effectively en- 
hances the robustness of a TTC. Fig. shows, for the 
full PPC-TTC model, time traces of the concentration of 
RpaA, the sum of the concentrations of the KaiC phos- 
phoforms that activate RpaA, and the sum of those that 
repress RpaA; for comparison, we also show the KaiC 
concentration in the TTC-only model. In the TTC-only 
model, the KaiC concentration varies slowly, and the con- 
centration threshold for kaiBC repression is crossed at a 
shallow angle. In contrast, in the PPC-TTC model the 
concentration of active RpaA switches rapidly between 
a value that is well above the activation threshold for 
kaiBC expression and one that is well below it. Such 
sharp threshold crossings minimize the effect of fluctua- 
tions in the concentration of the regulatory protein on the 
timing of the switch between kaiBC activation and kaiBC 
repression. RpaA's abrupt switching is a direct conse- 
quence of the sharp rise and fall in the concentrations of 
the KaiC phosphoforms that regulate RpaA. This in turn 




Time [h] Degradation Rate \i [1/h] 

FIG. 5: A PPC in combination with a TTC generates ro- 
bust rhythms over a wide range of degradation rates. (A) 
Time traces of the concentrations of active RpaA and of the 
KaiC phosphoforms that activate ("activator") and repress 
("repressor") RpaA for the PPC-TTC model (V = l^m 3 and 
fi = O.OShr" 1 ), and KaiC for the TTC-only model (V = l^m 3 
and fi = O.^hr" 1 ; for /i = O.OShr^ 1 the TTC does not os- 
cillate); (B) Correlation number of cycles as a function of 
degradation rate fj, for the PPC, TTC, and combined TTC- 
PPC model (V = l/im 3 ); for the TTC, n 1 /2 was computed 
for [C](£) rather than for p(t). The protein synthesis rates 
are varied with the degradation rates such that the average 
concentrations equal those used in the in vitro experiments 
(see Fig. E]). 



stems not so much from the oscillation in the total con- 
centration of KaiC, but from the phosphorylation cycle 
at the multiple phosphorylation sites of KaiC: The con- 
centrations of the various phosphoforms approach zero in 
certain parts of the cycle (Fig. [5^), even though the total 
KaiC concentration remains relatively large throughout 
(Fig. |4j\) . Indeed, the stability of the circadian rhythm in 
the PPC-TTC model is mostly determined by the KaiC 
phosphorylation cycle, which, as Fig. 03 shows, is highly 
robust. It thus seems that the phosphorylation cycle is 
the principal pacemaker, even though a transcription- 
translation cycle is needed to sustain it and can, in fact, 
raise its robustness above the limiting value at zero pro- 
tein synthesis and decay rate (Fig. [5j3) . 



III. DISCUSSION 

The evidence is accumulating that circadian rhythms 
are driven by both transcription-translation and protein- 
modification cycles not just in cyanobacteria but even in 
higher organisms [J, H§| . Our analysis suggests that both 
cycles are required to generate stable circadian rhythms 
in growing cells over a wide range of conditions. On the 
one hand, it is intrinsically difficult for a TTC to generate 
oscillations when the average protein decay time is long 
compared to the oscillation period. Indeed, the speed of 
protein degradation (and/or dilution) sets a simple up- 
per bound on the size of any oscillation; if only a small 
fraction of the proteins can be degraded in one oscillation 
cycle, then any modulation of the total protein concen- 
tration must necessarily have a low amplitude and so be 
very susceptible to noise. In cells whose doubling time is 
of order tens of hours (comparable to the 24 hour clock 
period) or longer, a TTC alone can thus function robustly 
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only at the expense of high rates of protein synthesis and 
active protein degradation. On the other hand, a pro- 
tein modification cycle must fail when protein turnover 
becomes faster than the clock period: Proteins are then 
typically degraded and replaced by new molecules in the 
wrong modification state before they can pass through 
the full cycle, destroying the oscillation. While a TTC is 
thus required to sustain a protein modification cycle at 
higher growth rates, a modification cycle can also greatly 
enhance the performance of a TTC, especially at lower 
protein synthesis and decay rates, where the modifica- 
tion cycle can be exploited to set the rhythm of gene 
expression. The protein phosphorylation cycle of KaiC 
is intrinsically very robust, because KaiC is present in 
relatively large copy numbers and because a KaiC hex- 
amer must go through many individual phosphorylation 
and dephosphorylation steps in each oscillation period. 
These two characteristics together are responsible for the 
reliable and pronounced oscillations in the concentrations 
of the various phosphoforms of KaiC (Fig. OV), which, 
in turn, drive the robust rhythm of kaiBC expression. 

Our combined PPC-TTC model is consistent with 
a number of experimental observations. It not only 
matches the average oscillations of the phosphorylation 
level and the total KaiC concentration in wildtype cells 
(Fig. UK), but also reproduces the observation that even 
in the presence of an excess of KaiA, the total KaiC con- 
centration still oscillates with a circadian period [l5[ (see 
SI). Yet, because quite a few elements of the TTC have 
not been characterized experimentally, our model of the 
TTC is necessarily rather simplified and phenomenolog- 
ical; not surprisingly, some observations thus cannot be 
reproduced. For example, our model predicts that the 
phase of the oscillation in KaiC abundance lags behind 
that of the phosphorylation level by a few hours, while ex- 
periments seem to show that these oscillations are more 
in phase 0, EH]. This may be due to our simplifying 
assumption that KaiC is produced as fully phosphory- 
lated hexamers. While phosphorylation of fresh KaiC has 
been reported to occur within 30 minutes Q, fragmen- 
tary evidence suggests that hexamerisation is slow [50l ]; 
one would then expect to detect KaiC monomers in ex- 
periments several hours before our simulations report the 
presence of fully functional hexamers. We believe, how- 
ever, that while our model is simplified, it does capture 
the essence of a coupled TTC and PPC: that is, it is 
built around a protein that not only regulates its own 
synthesis, but also undergoes a protein modification cy- 
cle, where the different modification states either activate 
or repress protein synthesis. Moreover, the basic trends 
we observe can largely be explained by the arguments 
of the preceding paragraph, which depend only on the 
comparison of a protein decay timescale with the oscilla- 
tion period; we thus expect that our qualitative results 
should apply to any system that exploits both a protein 



modification cycle and a transcription-translation cycle 
to generate circadian rhythms. 

In S. elongatus, the period of the circadian rhy thm is 
insensitive to changes in the growth rate [HI, l52j |. Since 
protein synthesis and decay rates, on the other hand, 
tend to vary with the growth rate [53| , the question arises 
whether the period of the oscillator as predicted by our 
model is robust to such variations. In the SI we show 
that a ten-fold increase in the synthesis and degradation 
rates of all proteins only decreases the period by 10-20%. 
The reason is that the period of the oscillator is mostly 
determined by the intrinsic period of the phosphoryla- 
tion cycle (Fig. [5^). The latter does not depend on the 
absolute protein synthesis and decay rates, although it 
does depend on the ratio of the concentrations of the Kai 
proteins [H, H, • Clearly, it would be of interest to in- 
vestigate how the ratio of the concentrations of the Kai 
proteins varies with the growh rate. 

Finally, how could our predictions be tested experi- 
mentally? The most important testable prediction of 
our analysis is that the PPC requires a TTC when the 
growth rate is high. Kondo and coworkers have demon- 
strated that a PPC can function in the absence of a TTC 
7]. These experiments were performed under constant 
dark (DD) conditions, which means that the growth rate 
was (vanishingly) small, protein synthesis was halted, 
and even protein decay was probably negligible, since 
the KaiC level was rather constant Q- These experi- 
ments are consistent with our predictions, but the criti- 
cal test would be to increase the growth rate while avoid- 
ing any circadian regulation of kaiBC expression. This 
means that the cyanobacteria must be grown under con- 
stant light (LL) conditions. Moreover, one would need to 
bring kaiBC expression under the control of a promoter 
that is constitutively active. However, most promoters 
of S. elongatus [33l . 13411 . and even many heterologous pro- 
moters from E. coli [35| . are influenced by the circadian 
clock. Nevertheless, a number of promoters have been 
reported that exhibit arhythmic activity 33j, and these 
might be possible candidates. A still more challenging 
experiment would be to express the phosphorylation cy- 
cle in the bacterium E. coli [5J]. Our analysis predicts 
that in growth-arrested cells, the phosphorylation cycle 
should be functional, while in normal growing E. coli 
cells, with cell doubling times of roughly 1 hour, the os- 
cillations should cease to exist. 
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I. THE MODEL 

A. Protein phosphorylation cycle 

The model for the KaiC protein phosphorylation cycle (PPC) is described by reactions (l)-(5) of the main text, 
which we repeat here: 

n ~ k f k 

C.^C,, C, + A <=* AC, -4 C i+ i + A, (SI) 

bi fe* b 
2fcf f _ kf f 

d + B <=* BC l5 BCi+B <± B 2 G l , (S2) 

kf b 2kf b 

xk^ f k^ { 

B x d + A ^ AB.C,, AB 2 Ci+A £± A 2 B 2 C 4 , (S3) 

kf b 2fc* b 

k k 

""pa , — a ps . 

Cj <=^ Gj+i, Ci Ci+x, (S4) 

fcd P= fcdps 

fcps fcps ^ 

B^Ci ^ B^Ci+X: A y B x Ci ^ A y B x Ci+i . (S5) 

Here, Ci denotes a KaiC hexamer in the active conformational state, in which the number i of phosphorylated 
monomers tends to increase, and Ci denotes a KaiC hexamer in the inactive conformational state in which i tends 
to decrease; A denotes a KaiA dimer, and B denotes a KaiB dimer. The reactions Ci C,; in Eq. [Si] model 
the conformational transitions between active and inactive KaiC; the second set of reactions in Eq. [HI] describe 
phosphorylation of active KaiC that is stimulated by KaiA; the reactions in Eqs. IS2I and [S3l model the sequestration 
of KaiA by inactive KaiC that is bound to KaiB; note that an inactive KaiC hexamer can bind up to two KaiA 
dimers; the reactions in Eqs. [S4] and [S5] model spontaneous phosphorylation and dephosphorylation of active and 



inactive KaiC. For a more detailed discussion of the model, we refer to 



ylat 



B. A PPC with constant protein synthesis and degradation rates 

In the main text, we also discuss the performance of a model that includes not only a PPC, but also synthesis and 
degradation of the Kai proteins. These reactions are given by 

h C 6 + 3B, h A, (S6) 

A,B,AC i ,B x C i ,A y B x C i 4 (S7) 

As explained in the main text, we assume that fresh KaiC is injected into the system as fully phosphorylated hexamers 
because phosphorylation of fresh KaiC proteins has been reported to be fast, i.e. occurring within 30 minutes 0. 
However, the precise choice for the phosphoform of fresh KaiC is not so important in this model; it does not affect 
the robustness of this model. Because the KaiB and KaiC proteins are both products of the kaiBC operon, we 
choose to model the production of both proteins as a single reaction. We note that while in the model with the 
transcription-translation cycle, discussed in section IIP! the delay in the synthesis reactions is critical, in the above 
model, where the Kai proteins are produced with rates that are constant in time, a delay would have no effect; the 
synthesis reactions are therefore modeled as simple, Poissonian birth reactions. 
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C. Deterministic PPC model with synthesis and degradation 



To verify that the disappearance of oscillations as the degradation rate is increased is not a purely stochastic effect, 
we consider the model of Eqs. (|S1|) — (|ST[) in the deterministic limit of infinite volume and protein number. In this limit, 
the concentrations of the different proteins evolve according to deterministic rate equations. We make two further 
simplifying assumptions: First, we replace the two-step binding of KaiB to with a trimolecular reaction that turns 
Ci directly into B 2 Ci, and making a similar change for binding of KaiA to the inactive branch. Second, we assume 
that binding and unbinding reactions are fast enough that they are effectively in steady state and thus explicitly keep 
track only of the concentrations of the various KaiC species; the concentrations of free KaiA and KaiB can then be 
inferred from conservation laws. The dynamical equations are then essentially the same as those given in Eqs. 44 - 47 
of our previous publication [l[ , with the addition of a linear decay term with rate [i for each species and of synthesis 
of C6 with rate /3 C : 



dt 

d\Ci 



dt 



d[B 2 C t 
dt 



e^JCwlT+^tC^lT-^r+^JICilT-af^QlT + of^G,] (S8) 
+AA,<» - f*[dh (S9) 
~k P s [Ci-i] + fcdps [Cm] - (fcp S + fcdps) [Ci] + af l [Ci] T - af h [C<] 

-hf([Bh 2 Ei[B 2 Q] T ) 2 [Q] + " f ~ K ' [B2 A C ; ]T - MCi] (S10) 

K, + [A] 2 

fcp S [B 2 C l -i]T + fcdps[B2C + i]T-(fcps + fcdps) [B 2 Ci]T 

+-? f ([B]x - 2EJB 2 Q] T ) 2 [C l ] - "V K ' [B ^ ]T ~ M[B 2 C l]T , (Sll) 

K, + [A] 2 



where the concentration of free KaiA, [A], is given by 

lA]+ ± p% +2 £ mcai-ni,... ( si2) 

Here [C,-]x is the total concentration of KaiC hexamers with i phosphorylated monomers in the active state, whether or 
not complexed with KaiA, i.e [Q]t = [Ci] + [AC,]; [B 2 Ci]x is defined similarly. The effective rate constants appearing 
in these equations depend on the concentration [A] of free KaiA and are defined in terms of the more microscopic rate 
constants as follows: The effective (de)phosphorylation rates on the active branch are af s = (fc ps Kj+fc p f[A])/(Kj + [A]) 
and af ps = Kik^ ps / (K$ + [A]). The effective flipping rates are given by erf ^ = /j-Ki/ (-fQ + [A]) and af 6 = b i: where fa 
and bi are the forward and backward flipping rates. The parameters kf 1 and kf h differ from kf l and kf b , respectively, 
in that the k's are rate constants for trimolecular reactions which are broken down into two successive bimolecular 
reactions in the stochastic simulations. The dissociation constants K, satisfy K^ = kf^/k^ 1 ; the K^ could be defined 
similarly in terms of forward and backward rates for KaiA binding to the inactive branch, but (just as with kf 1 
and kf b ) these rates would differ from those used in the stochastic simulations, so we choose instead to quote the 
dissociation constants directly. Following [lj, we choose values for the new parameters associated with trimolecular 
interactions such that time dependence of the phosphorylation ratio p(t) matches the average behavior of the stochastic 
model. 

To determine where oscillations disappear as /i is increased, we analyzed these equations using the XPPAUT 
implementation of the AUTO continuation package Q. We found that, for the parameter values given in Table [ 
the system undergoes a supercritical Hopf bifurcation at fj, = 0.0621 l /h, as noted in the main text. 



D. The PPC and TTC combined 



The full model of the main text consists of a PPC, a transcription-translation cycle (TTC), and a pathway that 
couples these two cycles. This model is described by the reactions of Eqs. IS1IIS5I for the PPC, together with the 
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Constant Value 
PPC (Eqs. 1-5 and lSHIS5)) : 
k pa ,k ?s 0.025 Vh 
kd ps ,k d ps 0.4 Vh 

fepf 1.0 l/h 

/i {10 -6 , 10" 5 , 10~ 4 , 10" 3 , 10" 2 , 10 _1 , 10} Vh 

bi 100 Vh 

fc^ f 1.72 ■ 10 10 i/M-h 

fc^ b {1,3,9,27,81,243,729} Vh 

fcf f 1.72 ■ 10 9 x {0.001, 0.1, 1, 1, 1, 1, 1} VM-h 

~kf b {10,l,l,l,l,l,l}i/h 

kf 1.72 ■ 10 9 x {0.01, 1000, 1000, 1000, 100, 0.001, 0.0001} i/M-h 

kt h {io, 1,1, 1,1,1, 10} Vh 

Deterministic PPC (Eqs. |S§1|ST2|) : 

Rf f 2.97 ■ 10 18 x {0.01, 1, 1, 1, 1, 1, 1} VM 2 h 

Rf h 100 x {10, 1,1, 1,1, 1,1} Vh 

Ki 3.37 • 10~ 25 x {oo, 100, 1, 1, 100, oo, oo} M 2 

RpaA activation (Eqs. 6 + IS13|) : 

K 8.6 • 10 9 i/M-h 

h 4.3 • 1 9 i/M-h 

TTC (Eqs. 7-9. IS14IS15|) : 

K 0.058 fiM 

A, (j,- 1 x 0.58 fM 

fit /i _1 x 0.29 /iM 

fi c from optimization for ([C]) = 0.58 /iM 

t 5h 
a T 0.5 h 

RpaA activation TTC-only (Eq. 11): 
k\ 1 VM-h 

fc? 100 VM h 

TABLE SI: The parameters used for the models of the main text. The degradation rate [i is a free parameter that we vary to 
explore different growth conditions. The numbers between the curly brackets in the right column correspond to the different 
KaiC phosphorylation states i in ascending order; values of co for indicate that a particular binding reaction is not allowed. 



following reactions for the TTC and the coupling between them: 

R + X^R + X, R + X^R + X, 

A[Ri m 4+[R]4) C 6 + 3B, 0^A, 

T±C7Y 

R, R, A, B, AC;, T^xCi, hyB x Ci A . 

Here, R and R denote the RpaA protein in its active and its inactive form, respectively. The X and X in Eq. (IS13[) 
denote any of the phosphoforms of KaiC that mediate the activation and repression of RpaA, respectively; in the 
next section, we discuss the dependence of the results on precisely which phosphoforms are chosen. The double arrow 
indicates a reaction with a Gaussian distributed delay with mean r and variance <j t . We thus assume that kaiBC 
expression is activated by RpaA, where the activity of RpaA is modulated by the PPC. In contrast, the expression of 
KaiA, KaiB, and RpaA is assumed to occur constitutively. 



(S13) 

^ R, (S14) 
(S15) 



E. Parameters 



Table [Si] gives the values of the kinetic parameters used in the stochastic simulations based on the kinetic Monte 
Carlo algorithm developed by Gillespie [4{. Unless otherwise noted, we choose the total concentrations of the Kai 
proteins to match common conditions for the in vitro reaction system 0, 0|: [A] T = 0.58 /iM, [B] T = 1.75 /iM, and 
[C] T = 0.58 ^M. 
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Model Activator 


Repressor 


Threshold K 


7l\/2 


a 


AC 2 ,AC 3 ,AC 4 , AC 5 


Aj / Bo;C2, • • 


.,A y B x C B 0.058 pM 


195 


b 


AC 3 ,AC 4 


A M B a; C2, . - 


.,A H B^C 6 0.058 /iM 


118 


c 


AC 3 ,AC 4 


AyBjjCa, . . 


.,A y B x C 4 0.058 


180 


el 


C3, AC3, C4, AC4 


AyB^Ca, . . 


.,A y B x G 4 0.029 /xM 


39 


e 


C X ,AC X , x £ {2,3,4,5} k y B x G 3 ,.. 


., A y B x G 4 0.029 (M 


48 



TABLE S2: Models with different output pathways from the PPC to the TTC. These models differ in the choice of phosphoforms 
that activate and repress RpaA, respectively. The maximal production rate /3 C has been modified such that the average concen- 
tration of KaiC is 0.58 fiM, as used in the in vitro experiments 5, 6J. In each case, m/2 is given for a volume V = 1 /im 3 and a 
decay rate fi — 0.03 h _1 . Model "a" is the full model from the main text. To make the simulations tractable, we neglected repres- 
sion of RpaA activation by KaiC phosphoforms that occur in negligible concentrations; consequently, the full list of phosphoforms 
that has the potential to repress RpaA is {B 2 C 6 , B 2 C 5 , AB2C5, A2B2C5, B 2 C 4 , AB2C4, A2B2C4, B 2 C 3 , AB 2 C 3 , A 2 B 2 C 3 , B 2 C 2 , 
AB 2 C 2 , A 2 B 2 C 2 }. Other parameters are given in Table [STl 



II. RESULTS ARE INDEPENDENT OF DETAILS OF THE OUTPUT PATHWAY 



In this section, we show that the precise choice of the KaiC phosphoforms that activate and repress RpaA is not 
critically important for the existence of oscillations. Tabic [S2] shows the different models that we have considered, 
and Fig. [ST] shows their time traces. It is seen that the time traces are very similar to those of the model in the main 
text, which is model "a" in Table [S2l The most significant difference can be observed for the time trace of RpaA in 
models "d" and "e". In these models, not only C^A activates RpaA, but also C^, thus KaiC that is not bound to 
KaiA. The concentration of C^A reaches zero during the dephosphorylation phase, and, as a result, the concentration 
of RpaA becomes zero during this phase in models "a"-"c". However, the concentration of Ca; does not reach zero 
during the dephosphorylation phase, and consequently, there is some residual activation of RpaA during this phase 
in models "d" and "e" . Nevertheless, RpaA activation during the dephosphorylation phase in these models does not 
manifest itself in the time traces of KaiC, because the concentration of active RpaA is still below the threshold for 
kaiBC expression. The oscillations of the phosphorylation level and total KaiC concentration are thus fairly similar 
in all models, although models "d" and "e" are less robust. 





FIG. SI: Time traces of the stochastic simulations of models with different output pathways from the PPC to the TTC (see 
Table [S2)l . The phosphorylation ratio (panel A) and the total concentration of KaiC (panel B) are fairly similar for all models 
studied. 
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FIG. S2: Contour plot of the correlation number of cycles of the PPC model of Rust et al. Q with production and degradation 
of Kai proteins with rates that are constant in time, as a function of the volume and the degradation rate. 



III. AN ALTERNATIVE PPC MODEL 



In the main text, we argue that the synergy between a transcription-translation cycle and a protein modification 
cycle is a generic feature of clocks that exploit both cycles. To support this claim, we have studied a model in which 
our model of the PPC is replaced by that of Rust et al. [7j. This model describes a phosphorylation cycle at the 
level of KaiC monomers, rather than KaiC hexamers as in our model. In the Rust model, each KaiC monomer 
cycles between an unphosphorylated state "U" , a singly phosphorylated state "T" where KaiC is phosphorylated at 
threonine 432, a doubly phosphorylated state "ST" where KaiC is phosphorylated at threonine 432 and serine 431, 
and a singly phosphorylated state "S" where KaiC is phosphorylated at serine 431 @, This cycle is described by 
the reactions 

U O T, T o ST, ST o- S, S o- U (S16) 

with reaction rates given by Eq. 5 of the supplementary material of Rust et al. [7]. These rates depend on the 
concentration of free KaiA, which is sequestered by KaiC in the S-state. We model KaiA sequestration explicitly: 

S + AoAS, AS + Ao A 2 S. (S17) 

We picture sequestration to be fast and we picked a forward rate of 0.58 Y^M-s and a backward rate of 1 /s for both 
equations above. Dephosphorylation of KaiC in the S-state might occur even when KaiA is bound, in which case the 
KaiA protein is released from the complex. We define the output signal as 

(f) [T] + [ST] + [S] + [AS] + [A 2 S] 
P[ ' [U] + [T] + [ST] + [S] + [AS] + [A 2 S] ' 1 ' 

which resembles the phosphorylation ratio in the case where we cannot distinguish between singly and doubly phos- 
phorylated KaiC. The denominator in the above expression is also the total KaiC monomer concentration. We use the 
same concentrations as Rust et al., [KaiA] = 1.3 (active KaiA monomers) and [KaiC] = 3.4 /zM (KaiC monomers), 
and simulate this model using the Gillespie algorithm 

When we simulate this model for a volume V = 1 /xm 3 , we find a period of L = 21.4 h and a decay constant for the 
autocorrelation function of = 128 h. The corresponding correlation number of cycles is ni/2 = 30, which is lower 
than that observed experimentally, n l / 2 — 166 ±100 @, and lower than that of the PPC developed by us [l|, for which 
rc.i/2 ~ 200. This is because the model of Rust et al. features a phosphorylation cycle at the level of KaiC monomers 
rather than KaiC hexamers as in our model. The concomitant reduction in the total number of phosphorylation steps 
in the cycle reduces the robustness in the model of Rust et al. Q ■ 



A. The Rust model with constant protein synthesis and degradation 



To study the behavior of the model of Rust et al. [7] under conditions in which cells grow and divide, we have to 
include protein degradation and make up for this by protein synthesis. As in the main text, when we vary the protein 
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FIG. S3: Time trace of the Rust model [3] extended to include a TTC for V = 1 /im 3 . Both the total concentration of KaiC 
and the phosphorylation level of KaiC show stable oscillations. The messenger protein RpaA concentration is shown tenfold 
and crosses the threshold of K — 0.174 /iM (w 100 molecules) reliably. 



degradation rates, we adjust the protein synthesis rates such that the average protein concentrations are unchanged 
and similar to those used in the in vitro experiments [H, Q . 

Figure [S2l shows the results for this model. They are qualitatively the same as those of Figure 3 of the main text: 
the robustness decreases with decreasing volume and increasing degradation rate. Hence, not only in our model but 
also in that of Rust et al. protein degradation can cause the oscillations to disappear. This supports our claim 
that a protein modification oscillator cannot function on its own when the cell's growth rate is high enough. 



B. The Rust model with a TTC 

We will now show that a TTC can resurrect the PPC of Rust et al. [7J. We model the TTC as: 

*MVJ£ + M-) gT) % A 

r,a r 

R + T^R + T, R + A X S h R + A^S, 

A,R,R,U,S ) A a; S,T ) ST4 0. 

where the first line describes the production of proteins to counteract their degradation and the second line summarizes 
the RpaA signaling pathway, where KaiC that is phosphorylated at the T site activates RpaA and KaiC that is 
phosphorylated at the S site represses RpaA activation. The parameters in this model are f3 c = 1.13 ** M /h, K = 
0.058 fM, /3 a = 0.13 mM/ Sj (3 r = 0.058 mM/ 8) k a = k, = 0.1 i/W-s, and the delay is r = (3 ± 0.3) h. Figure El shows the 
results of this model at a volume V = 1 /im 3 . Analysing the autocorrelation function, we find a period L = 22.5 h and 
a correlation decay time of t = 1832 h leading to a correlation number of cycles of ni/ 2 = 402. Clearly, a TTC can 
also resurrect the PPC of Rust et al. Q, supporting our claim that the qualitative results of the main text should 
apply to any biological system that exploits both a protein modification cycle and a protein synthesis cycle to generate 
circadian rhythms. 



^ R, (S19) 

x 6 {0,1, 2}, (S20) 
(S21) 



IV. THE EFFECT OF BURSTS 



In the model of the main text, described in section U of this SI, we have concatenated transcription and translation 
into one gene-expression step. Moreover, we have ignored promoter-state fluctuations. Allowing for the explicit 
formation and translation of mRNA [10], as well as for slow promoter-state fluctuations [lj [l2|], could lead to bursts 
in protein synthesis, which are expected to lower the robustness of the TTC. This could potentially lower the stability 
of the clock. To address this, we have performed simulations of a model in which KaiB and KaiC are produced in 
bursts. We assume that 5 KaiC hexamers and 15 KaiB dimers are formed in each gene expression reaction (rather 
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FIG. S4: The correlation number of cycles as a function of the degradation rate fx and the volume V for a TTC-PPC model 
that exhibits bursts in gene expression; upon a gene expression event, 5 KaiC molecules are produced and 15 KaiB dimers. 



than the 1 and 3 of Eq. IS 14[) : this corresponds to typical burst sizes observed experimentally [10] in E. coli. Eq. IS 141 
is thus replaced by: 

WRi m 4+[R]4) 5C 6 + 15B. (S22) 

T±ov 

Fig. [S4] shows the resulting phase diagram. As expected, it is similar to Fig. 4 of the main text, which shows the 
results of the PPC-TTC model without bursts in gene expression. The robustness of the model with bursts is lower, 
but not very much so: ni/a = 150 for the model with bursts versus n\i% = 195 for the model without bursts shown 
in the main text (// = 0.03hr _1 and V = 1/xm 3 in both cases). We believe that this relatively small reduction in the 
clock's stability is due to the stabilizing effect of the PPC. 



V. THE FULL MODEL WITH VOLUME GROWTH AND BINOMIAL PARTITIONING 

Living cells constantly grow and divide, and proteins thus have to be synthesized to balance dilution. In the 
main text, we argued that the principal effect of dilution is to introduce an effective degradation rate set by the 
cell doubling time. Here, we show that this is indeed the case: we study a model in which growth, cell division 
and binomial partitioning of the proteins upon cell division are modeled explicitly and show that its qualitative 
behavior is similar to the model of the main text, in which the volume is held constant and protein degradation occurs 
at a rate that is constant in time. 

The model we consider here is the combined PPC-TTC model presented in the main text, but with the degradation 
reaction, Eq. 9, replaced by a scheme in which the bacterial volume V grows exponentially as 

V(t) = V e t! %t, (S23) 

where T$ denotes the doubling time after which the volume reaches twice its minimum V Q and cell division is triggered. 
Division includes dividing the volume by two, partitioning the proteins binomially [1J|, and deleting events on the 
queue of the delay associated with the KaiC production reaction with a probability of 0.5 for each daughter cell. 
To compare the results of this model with those from the main text, we take = ln2//i, where fi is the protein 
degradation rate of the model in the main text; if proteins were to decay only by dilution in a cell with a doubling time 
Td, then /j. would be the effective protein degradation rate; if proteins are also degraded actively, then /i = In2/Td is 
a lower bound on the actual degradation rate. 

Fig-ESshows time traces of the total KaiC concentration and the KaiC phosphorylation level for this refined model. 
It is seen that the oscillations of the total KaiC concentration are more noisy than those in the model in which the Kai 
proteins are produced and degraded with rates that are constant in time (Fig. 4A of the main text). Clearly, binomial 
partitioning is a major source of noise, with the random removal of items from the queue associated with the KaiC 
synthesis reaction being the largest source of noise. Nonetheless, the oscillations of the KaiC phosphorylation level 
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FIG. S5: Time traces of the full model, given by Eqs. 1-9 of the main text, modified to take into account cell division and 
binomial partitioning. Here, the cell divides when the volume reaches V m = 1.38 /im 3 , which occurs every 20 h, as indicated 
by the arrows above the graph; a cell doubling time of 20 h corresponds to an effective degradation rate of p = 0.035 Yh. The 
average volume is V = 1 /jm 3 . (A) Time traces of the phosphorylation level and the total KaiC concentration. (B) Time traces 
of the total KaiC concentration, the KaiC copy number and the volume. Please note that the time trace of the total KaiC 
concentration is hardly affected when cell division happens to occur during the degradation phase, while it has a relatively 
large effect when cell division happens to occur during the KaiC production phase; this is because the removal of items from 
the queue associated with the KaiC synthesis reaction effectively lowers the synthesis rate; indeed this explains the change in 
slope in the KaiC concentration during the production phase. 



are much less affected, with the correlation number of cycles being n 1 / 2 = 88. Indeed, while this model combining 
a TTC with a PPC is fairly robust, an oscillator with exponential volume growth and binomial partitioning built 
upon a TTC alone, is not stable. This supports our statement in the main text that a PPC can strongly enhance the 
robustness of a TTC. In a future publication, we will systematically study the effect of bursts in gene expression and 
binomial partitioning. 



VI. ALTERNATIVE TTC MODELS 



We not only studied the TTC model discussed in the main text, but also the simplest possible TTC model, namely 
one in which KaiC directly represses its own synthesis (with a delay): 

^lK +[C]i) C, C A 0. (S24) 
T;oy 

The result of this model is shown in Fig. [S6j It is seen that this TTC model behaves very similarly to that of the main 
text (compare Fig. 5B): it shows robust oscillations only for high degradation rates. Also including enzyme-mediated, 
zero-order protein degradation [l4j does not qualitatively change this behavior. This supports our argument that the 
TTC has an intrinsic difficulty in generating oscillations when the doubling time is long and there is no strong active 
degradation. 

In the simulations of Fig. 5B of the main text and Fig. [56] of the SI we vary the protein synthesis rate /3 C and 
the protein-production delay r such that the average KaiC concentration remains constant at the in vitro level [f| 
and the period remains constant at 24 hours. This procedure means that as the degradation rate is increased, not 
only the synthesis rate has to be increased, but also the delay. One may wonder what would happen if the delay 
were kept constant and only the synthesis rate was adjusted, such that only the period remains constant at 24 hours; 
the average KaiC concentration is thus allowed to change. To address this, we show in Fig. [57] the results of a 
deterministic calculation based on the mean-field, chemical rate equations. Panel A shows as a function of the protein 
degradation rate /z, the KaiC synthesis rate f3 that keeps the oscillation period at 24 hours, for different values of the 
delay. It is seen that if the delay would be kept constant at 5-6 hours, as in the full model, the synthesis rate would 
have to be increased dramatically when the degradation rate is increased, to keep the period at 24 hours. This would 
raise not only the average KaiC concentration, but, more importantly, also the protein synthesis rate] this would 
increase the energetic cost of sustaining a TTC dramatically, as shown in Fig. IS7B . In the main text, we keep not 
only the clock period constant, but also the average KaiC concentration. This means that when the degradation rate 
is increased, not only the synthesis rate, but also the delay is increased; however, even in this scenario, the cost of 
synthesizing proteins still increases markedly (Fig. IS7B ). 
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FIG. S6: A comparison of the performance of the PPC-only model (Eqs. I S 1 1 — IS7 [ l. the full model including both PPC and 
TTC (Eqs. ED -EH and Eqs. IS131 - [ST5]) . and the simplest possible TTC model, namely one in which KaiC directly represses 
its own synthesis (see Eq. IS24|) . Compare to Fig. 5B of the main text. 



VII. PERIOD AS A FUNCTION OF CELL VOLUME AND PROTEIN DEGRADATION RATE 

Fig. EH shows the period of the oscillation of the KaiC phosphorylation level in the full model of the main text 
(Eqs. 1-9) as a function of the cell volume and the protein degradation rate. As before, the protein synthesis rates are 
adjusted such that the average protein concentrations are constant and similar to those used in the in vitro experiments 
5, 6]. It is seen that the dependence of the oscillation period on the cell volume and protein degradation rate is rather 
weak. This is because the rhythm of the clock is dictated by the PPC, which is insensitive to the absolute rates of 
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FIG. S7: The effect of changing the delay and the protein synthesis rate in a simple TTC based on Eq. IS24I (A) The KaiC 
production rate /3 that keeps the oscillation period at 24 hours, as a function of the KaiC degradation rate, for different values 
of the delay r. (B) The average synthesis rate (i.e. the average number of KaiC molecules produced per clock period), which is 
a measure for the cost associated with protein production during one period, as a function of the degradation rate. The model 
is a simple TTC in which KaiC directly represses its own synthesis, given by Eq. IS24I The calculations were performed using 
the deterministic, mean-field chemical rate equations, because performing stochastic simulations with high production rates is 
prohibitively expensive. The crosses denote the points where not only the oscillation period is 24 hours, but also the average 
KaiC concentration equals that used in the main text, which is the KaiC concentration used in the in vitro experiments [1, 0- 
The lines end for low degradation rates because the oscillations cease to exist. Panel A shows that to keep the clock period at 
24 hours, the synthesis rate has to increase with the degradation rate, if the delay is kept constant. This not only raises the 
average KaiC concentration, but also increases the energetic cost of sustaining the TTC, as shown in panel B. 
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FIG. S8: Period of the oscillation in the KaiC phosphorylation level in the full model of the main text as a function of cell 
volume and protein degradation rate. When the degradation rate is varied, the protein synthesis rates are adjusted such that 
the average protein concentrations are constant and similar to those used in the in vitro experiments 0, 0|. The period is 
essentially independent of the volume, and exhibits only a weak dependence on the degradation rate /i. 



protein synthesis and decay. We note here that the period of the oscillation, as well as its amplitude, would change 
if the ratio of the concentrations of the Kai proteins were changed. While the dependence of both the amplitude and 
theperiod of the in vitro PPC on the ratio of the concentrations of the Kai proteins has been characterized in detail 
0; HI j the dependence of the in vivo oscillator on their ratio has not been studied experimentally. 



VIII. RHYTHMS OF KAIBC EXPRESSION WHEN KAIA IS OVEREXPRESSED 



Kitayama et al. have shown that kaiBC expression oscillates with a circadian period in the presence of an excess of 
KaiA, although it is not clear whether these oscillations are sustained or damped [Hj]. Our PPC-TTC model of the 
main text, which is described by Eqs. IS H IS5I and Eqs. IS 131 — [STol above, generates oscillations in kaiBC expression 
with a period of 24 hours when kaiA is overexpressed threefold, as shown in Fig. IS9A . This figure shows that the 
phosphorylation level also exhibits weak oscillations, which are not seen experimentally; this is due to the fact that 
our PPC model neglects phosphorylation of inactive KaiC, which is known to occur at high KaiA concentrations 
To rectify this, we have extended our model to include KaiA-stimulatcd phosphorylation of inactive KaiC, using the 
same reactions as those used for active KaiC; specifically, we add to the reactions of Eqs. [STHS5land Eqs. IS13I — [S~15l 
the following reactions: 

A.yB.C, + A <± A^C.A 4 A y B x C l+1 + A, (S25) 

for each phosphorylation level i; the rate constants equal those of the corresponding reactions of active KaiC, except 
that the KaiA-KaiC association rate is reduced by a factor 100. All other rate constants are as in Table [Si] We also 
include auto-activation of RpaA via 

R % R, (S26) 

with /c™ = 25hr _1 . Auto-activation of RpaA becomes necessary because the concentrations of the KaiC phosphoforms 
that activate RpaA become very low when KaiA is in excess. (The freshly injected KaiC hexamers don't make it 
to the bottom of the phosphorylation cycle because of the excess KaiA.) This model not only matches the in vitro 
observation that, when an excess of KaiA is added during the dephosphorylation phase, the phosphorylation level of 
KaiC rises immediately Q, but also reproduces the in vivo oscillations of the total amount of KaiC when KaiA is 
overexpressed [15], as shown in Fig. IS9B. 



(A 



(B) 



FIG. S9: The PPC-TTC model predicts oscillations in kaiBC expression when kaiA is overexpressed. (A) The original PPC- 
TTC model of Eqs. ISll-[S5land Eqs. IST3l-fST5l(B) The modified PPC-TTC model, given by the original model of Eqs. [SlHS5l 
and Eqs. IS 131 — [S15l but extended to include KaiA-stimulated phosphorylation of inactive KaiC (Eq. IS25I) and auto-activation 
of RpaA (Eq. IS26I) . Note that the total amount of KaiC exhibits pronounced, albeit somewhat arhythmic, oscillations, while 
the phosphorylation level is high and fairly constant. For both panels, the protein degradation rate is /i = O.lhr" 1 and the 
volume is V = 2.97/^m 3 ; kaiA is overexpressed by a factor of 3. 



IX. MEASURING THE ROBUSTNESS 



In this section, we discuss how we calculate the correlation number of cycles rii/2 for our various models. We begin 
with some theoretical background: Consider a phase variable ip(t) that increases with an average frequency u) and is 
also subject to noise. Its time evolution can be written as 

M^=a; + £(t) with (£{t)£(t')) =a 2 5(t-t') . (S27) 

Here is Gaussian white noise of strength er 2 , and (o) indicates averaging over different realizations of the noise. 
Integrating the equation with the initial condition <p(0) — yields 

ip(t) =ujt + W(i) , (S28) 

where W(t) is a Gaussian random variable with mean zero and variance 

((<p(t) -Lot) 2 ) = cr 2 t. (S29) 
From this, we can construct an oscillating signal 

x{t) — x + a ■ sinip(t) (S30) 
with mean x and amplitude a. The autocorrelation function of (|S30I) is then 

C(tO = ^ff^=e-^KI.cosK), (S31) 

(6x(ty) 

where Sx(t) = x(t) — xq is the deviation of the signal from its average value xo- We thus expect that the correlation 
function is a sinusoid modulated by a single exponential decay. 



A. Incorporating Amplitude Noise 

Because the phase is the only "soft" direction of the dynamical system, one generically expects any noisy oscillator 
at long times to act similarly to the simple model of a phase oscillator just described. Nonetheless, a more realistic 
description would also include a fluctuating amplitude. We incorporate this behavior phenomenologically by including 
a time dependence in the amplitude: 



a -)• a(t) = a Q + £(„)(*) 



x(t) = xo + (a + £(„)(*)) • sin(wi + W(t)) 



(S32) 
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where £Ww denotes a Gaussian white noise process of strength a 2 and therefore neglects correlations in the amplitude 
fluctuations. We can again calculate the correlation function analytically. By definition, we have C(0) = 1 for t = 0, 
but because a(t) now contains a white noise term, C(t) jumps discontinuously to a smaller value for any t > 0. One 
finds 

w- ^tm™ -^'-™-"™- ,,>0 ' (S33) 

Considering the more natural case of a finite correlation time in the amplitude noise, the picture will only change 
slightly: instead of jumping from 1 to v at r = 0, the envelope will undergo a smooth transition for this envelope 
involving two time scales: a short time scale of order the correlation time of the amplitude fluctuations and a much 
longer time scale associated with the phase diffusion. In practice, we found that including amplitude fluctuations did 
not significantly change our estimate for n\j%. 



B. Computing the Correlation Number of Cycles 



To calculate the correlation number of cycles rii/2, we begin by using our simulation results to estimate the corre- 
lation function C(t). After an initial equilibration phase of 500 hr, we do simulations for 50,000 hr. From these, we 
extract N = 500, 000 values Xi for the time trace x(t) at equidistant points in time, which we can then use to calculate 
the correlation function at times separated by an interval At = 0.1 h: 

^ N-i 

Xi = x(t + i ■ At), C(i ■ At) = ^ N _.y ^ 6x i 5x i+i ■ ( S34 ) 

Here x(t) is either the phosphorylation level p(t) or the total concentration of KaiC hexamers. p(t) is used for all 
cases except for models where there is only a TTC and a phosphorylation level cannot be defined. 

Once C{t) has been computed, it can be fitted to the form (|S33|) to determine the free parameters v, a 2 , and w. In 
practice we perform the fit using gnuplot, which in turn implements a Marquardt-Levenberg algorithm. 

Finally, it remains to translate the fitted parameter values into an estimate of ni/ 2 . Eguchi et al. define ni/ 2 as 
the number of cycles after which the standard deviation of the phase is n [l6j]. Thus, we have 

yj ([<p(L ■ n 1/2 ) - art]*) = tt => n 1/2 = ^, (S35) 

where L = 2tt/oj denotes the period and we have used Eq. (|S29|) to solve the left equation for n 1 / 2 . The value of n X j 2 
is then obtained by substituting the fitted values for a 2 and to. Using this method we can reliably measure n.1/2 from 
1 to 1,000. The upper bound is given by fact that we only compute the correlation function C(t) up to t — 3, 000 h 
and therefore correlation functions which decay on much longer timescales are difficult to detect. 

Concerning the error bar on the computed ni/ 2l it should be noted that there are two distinct sources of error: one 
due to statistical fluctuations, and one due to systematic errors, e.g. that the correlation function cannot be fitted to 
Eq. IS33I To estimate the former, for certain parameter values we repeated the procedure described above 32 times, 
i.e. we performed 32 independent simulations and computed the mean and the standard deviation of the set of 32 
independently estimated values for n 1 / 2 . For V = 1/im 3 and fi = 0.025hr _1 , this gave ni/ 2 = 168 ± 31 (mean ± std. 
deviation) for the combined PPC-TTC model. To address the second type of error, we also computed n 1 / 2 using two 
different methods. One is the method of Eguchi et al., which is based on examining the times at which the oscillation 
reaches its maximum in each cycle, and thus does not assume any particular form like Eq. (|S30[) for the oscillating 
signal 16]. The other method involves computing the width of the dominant peak in the power spectrum of the 
time traces. The three methods gave similar error bars and values for n\j 2 that agree within the error bar. While 
all methods gave the same result, we found the method based on the correlation function (Eq. IS33|) more robust for 
noisy oscillations at low cell volumes. 
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